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ABSTRACT 

I , 

A computer  model  Is  developed  for  simulation  of  two  acoustic 
signal  processors:  the  real  time  processor  (RTP)  and  the  frequency  ^ 

domain  processor  (FFT).  | j 

i 

Simulations  are  performed  to  assess  the  relative  merits  of  the 
two  processor  techniques  and  comparisons  of  the  possible  performance  i 

of  the  two  systems  are  presented.  | 
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COMPARISON  OF  REAL  TIME  AND  FFT  PROCESSING 


ALGORITHMS  FOR  ACOUSTIC  DATA 

1.0  SUMMARY 

Use  of  acoustic  signals  emanating  from  possible  targets  as  a 
source  for  detection,  location  and  tracking  to  these  targets  is  of 
Interest  to  the  army  RPV  program. 

Acoustic  signals  generated  by  the  type  of  targets  usually  sought 

I 

generally  reside  in  the  low  audio  frequency  range.  The  ambient  noise 
background  is  often  considerable  in  this  frequency  range  and  it  is 
necessary  to  provide  a somewhat  sophisticated  algorithm  for  processing 
the  acoustic  signals. 

Algorithms  for  processing  these  signals  and  determining  the 
relative  bearing  of  the  targets  have  been  developed  through  previous 
efforts  at  other  agencies.  These  algorithms  generally  may  be  classi- 
fied as  one  of  two  basic  methods:  Real  Time  processing  and  Frequency 
Spectrum  Processing  (generally  a Fast  Fourier  Transform  (FFT)  applica- 
tion) . 

Although  the  FFT  processing  techniques  have  held  the  advantage 
until  recently  in  so  far  as  accuracy  is  concerned,  recent  state-of- 
the  advances  in  software  processing  and  Integrated  circuit  chips 
(IC's)  lead  one  to  believe  that  a Real  Time  processing  (RTP)  algorithm 
for  accuracy  and  at  the  same  time  update  data  faster  than  the  FFT 
algorithm.  Furthermore  the  RTP  algorithm  generally  requires  conside- 
rably less  hardware  than  the  FFT  algorithms. 
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This  report  documents  an  Investigative  study  of  the  two  basic 
algorithms,  RTF  and  FFT,  and  presents  the  results  of  a computer 
simulation  study  of  their  relative  performance. 

Under  some  circumstances  the  RTF  is  far  superior;  for  instance 
in  the  situation  where  the  signal  frequency  varies,  such  as  an  engine 
changing  speeds.  Both  algorithms  suffer  inaccuracies  for  target  angles 
off  the  end  of  the  array  (small  angles)  but  this  presents  no  real 
problems  for  a tracking  application  since  the  target  will  be  kept  in 
a broadside  situation  to  the  array  (large  angles)  by  the  tracker. 

In  general  it  is  concluded  that  a RTF  algorithm  is  a very  meri- 
table  approach  and  should  be  realized  in  a hardware  prototype  form  to 
use  in  comparitive  studies  with  an  FFT  algorithm  hardware  prototype. 

This  research  work  funded  a graduate  student  working  towards  an 
MSEE  degree  and  a paper  describing  the  results  of  this  research  is 
being  prepared  for  publication  and  presentation. 

Much  appreciation  is  expressed  to  Captain  H.  Norckauer,  Richard 
Currie  and  Steve  Golden  of  the  Advanced  Sensors  Directorate  (IR  Branch) 
at  Redstone  Arsenal  for  their  interest  and  suggestions. 
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2.0  INTRODUCTION 

Acoustic  signals  emanating  from  targets  such  as  tanks  have  been 
studied  and  their  characteristics  noted  for  use  In  detecting  and 
tracking  devices.  One  characteristic  assumed  for  far  field  distances 
Is  the  shape  of  the  acoustic  wavefront  as  It  propagates  outward.  The 
wavefront  Is  assumed  to  be  planer,  that  Is  It  Is  assumed  to  be  an 
Infinite  plane  wave  (or  a spherical  wavefront  of  Infinite  radius). 

Thus  the  normal  vector  to  the  wavefront  points  to  the  source  of  the 
wavefront,  which  Is  assumed  to  be  the  target  location.  This  normal 
vector  describes  an  angle  with  respect  to  the  axis  of  a two  microphone 
array  which  may  be  determined  by  measuring  the  phase  difference  of  the 
waveform  as  It  crosses  the  two  microphones. 

The  angle  desired  Is  6,  Illustrated  In  figure  1.  This  angle  can 
be  derived  by  several  methods.  If  either  real  time  waveforms  or  FFT 
algorithms  are  utilized  and  a periodic  source  assumed,  then  the  angle 
Is  found  by  noting  that 


where:  d ■ the  array  spacing 

and  h - T 

with  Vgj  the  velocity  of  sound  and  t the  time  difference  In  arrival 
of  wavefront  at  microphones  1 and  2 and 

tCl/freq.) 

^ (360^)  • 
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Wavefront  travels  at  velocity  > llOOFT/SEC 

0 • Microphones 

t “ time  to  travel  distance  h at  1100  ft/sec  - 
(p  Shift  from  Mike  1 to  Mike  2 - 360“  x ^ - - 


T IIOO(T) 

T “ Period  of  the  frequency  components  of  acoustic  wavefront 
(assumed  to  be  - 100  Hz  T - .01  sec.) 


e ■ Target  Angle 


cos 


j 


-1  i|i  Shift  X 1100  X T 


4 X 360 


■] 


Example:  If  p Shift  Is  calculated  to  be  100.28218“  for  100  Hertz 
component  with  d “ 4 feet  for  mike  spacing.  Then; 


0 - 


-1 

cos 

*1 

100.28218  X 1100  X .01 

* 

4 X 360 

40.000“ 


Figure  1.  Two  Microphone  Acoustic  Array 


where  H Is  the  phase  difference  of  the  waveforms  at  1 and  2,  and 
freq.  Is  the  frequency  of  the  waveform.  Thus 
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-1 

cos 


1I(l/freq.)V 

d(W^ 


Inaccuracies  arise  from  the  following  assumptions: 

A.  The  acoustic  wavefront  is  assumed  to  be  a plane  wave  over  a 
distance  of  a wavelength  and  to  have  a normal  vector  which  passes 
through  the  centroid  of  the  sound  source  and  through  the  center  of 
the  microphone  pair. 

B.  The  acoustic  wavefront  Is  assumed  to  be  a periodic  single 
frequency  wavefront.  (Note  the  use  of  the  frequency  assumed  In  the 
equation  for  6 . ) 

C.  The  velocity  of  sound  Is  assumed  to  be  constant  and  known  so 
that  the  distance,  d,  represents  a certain  percentage  of  the  wave- 
length. 

Items  of  Note: 

(1)  The  spacing  of  the  array  Is  assumed  to  be  such  that  no  more 
than  180°  phase  shift  will  occur  when  the  wavefront  passes  the  array. 
Otherwise  the  array  cannot  differentiate  which  hemisphere  the  sound 
originates  from  with  respect  to  a line  drawn  perpendicular  to  the 
line  containing  the  two  microphones.  Thus  the  array  spacing  Is 
determined  In  part  by  the  highest  frequency  one  expects  to  detect. 

(2)  A second  array  Is  necessary  to  determine  which  quadrant  the 
sound  originates  from  since  a single  two-mlcrophone  array  does  not 
determine  the  Quadrant. 
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It  is  feasible  to  determine  the  velocity  of  sound  across  the  array 
by  measuring  the  actual  temperature  as  the  phase  Is  being  measured. 

Item  C.  This  point.  Incidentally,  Is  a problem  for  both  real  time 
systems  and  for  FFT  systems.  Measurement  of  the  temperature  at  the 
array  center  provides  sufficient  correction  for  the  calculation  of  the 
angle  and  It  is  noteworthy  to  mention  that  It  Is  not  necessary  to 
determine  the  temperature  along  the  path  between  the  source  and  the 
detector  since  the  velocity  of  the  waveform  is  incidental  to  the  calcu- 
lations except  during  the  time  of  travel  from  the  first  microphone  In 
the  array  to  the  last  microphone  in  the  array.  The  assumption  that 
the  velocity  of  the  wavefront  Is  constant  is  valid  since  the  array 
size  Is  small  and  we  are  Interested  again  In  the  velocity  only  as  the 
wavefront  crosses  the  array. 

Both  the  real  time  and  the  FFT  processors  will  suffer  inaccuracies 
due  to  non-periodic  components  of  the  input.  Item  B.  In  the  FFT  pro- 
cessor, however,  the  noise  and  the  varying  spectral  content  of  the 
signal  are  treated  as  periodic  signals  over  each  sampling  period.  As 
a result,  these  variations  are  averaged  out  to  some  extent  and  the 
spectral  content  of  the  target  is  enhanced  with  respect  to  the  noise 
background . 

The  real  time  processor  does  not  inherently  average  over  any 
period  and  as  a result  the  continuous  output  of  phase  difference  infor- 
mation varies  considerably  due  to  noise  and  spectral  variation  of  the 
signal.  In  a real  time  processor  as  designed  by  this  author,  however, 
the  filter  bandwidth  may  be  considerably  smaller  than  that  of  the  FFT 
processor  which  does  tend  to  reduce  the  noise  effects  considerably. 
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It  Is  also  beneficial  to  devise  some  sort  of  averaKln^  scheme  to  use 
with  the  real  time  processor  to  Improve  the  accuracy  and  to  further 
reduce  the  effects  of  noise. 

The  follovlnf;  describes  a software  approach  (now  feasible  with 
LSI  and  microprocessor  technology)  which  provides  an  "averaging"  effect 
for  a real  time  process  while  allowing  faster  processing  time  than  an 
FFT  process. 

The  frequencies  of  Interest  are  typically  In  the  range  of  100 
Hertz.  Using  real  time  waveforms  and  a phase  locked  loop  phase  com- 
parator the  phase  angle  may  be  measured  directly.  It  Is  proposed  to 
sample  the  phase  at  certain  Intervals,  say,  .05-second  Intervals  for 
ten  samples.  Let  each  phase  sample  be  compared  to  the  previous  and  If 
It  agrees  with  another  sample  within,  say,  ten  degrees,  average  the 
two  and  store.  If  It  does  not  agree,  store  for  further  comparison  as 
the  possible  target  but  with  low  priority.  In  this  manner,  the  ten 
sample  processor  will  discard  large  errors  and  average  the  small  errors. 
After  ten  samples  are  used,  an  updating  may  be  accomplished  by  simply 
keeping  a running  average  and  discarding  the  earliest  data  sample. 

While  It  Is  true  that  this  procedure  will  tend  to  average  out  valid 
target  angle  changes  In  .5  second  Intervals,  the  error  at  one  mile  would 
be  1.09°  at  the  worst  case  and  for  shorter  ranges  should  decrease  since 
the  bearing  change  would  be  kept  relatively  small  by  the  tracking  system. 

A second  and  major  problem  arises  when  the  target  Is  an 
engine  which  Is  changing  speed.  This  generates  a signal  which  Is 
changing  In  frequency.  The  FFT  processor  which  uses  long  time  averaging 
(from  0.5  second  to  2.0  second  intervals)  will  detect  this  type  of 
signal  as  energy  existing  over  a wide  frequency  spectrum  and  will 
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exhibit  large  errors  In  the  angular  calculations  of  such  a target. 

On  the  other  hand,  the  real  time  processor  uses  very  short 
looking  Intervals  and  does  not  see  appreciable  frequency  change.  This 
results  In  fairly  accurate  angular  calculations  especially  for  angles 
around  the  null  zone  (6  ■ 90“). 

Results  are  presented  which  Illustrate  \’ery  nicely  the  ability 
the  real  time  processor  to  handle  this  situation  and  the  Inability 
of  the  FFT  processor  to  deal  with  this  problem. 

One  possible  solution  for  the  FFT  processor  to  deal  with  this  Is 
to  use  a frequency  tracking  filter  but  this  significantly  adds  to  the 
cost  and  complexity  of  the  system  and  has  not  been  Investigated  to 
determine  any  benefits  In  performance. 

Assumption  that  the  normal  vector  to  the  acoustic  wavefront  passes 
through  the  centroid  of  the  sound  source.  Item  A,  Is  not  necessarily  a 
valid  assumption.  Ripples  In  the  wavefront,  tilt  due  to  propagation 
effects,  multipath  reception  all  combine  to  create  errors  In  the  normal 
vector  pointing  direction.  These  errors,  however,  are  typically  small. 
In  the  range  of  5*  or  less,  and  are  not  troublesome  during  the  tracking 
mode  until  one  reaches  a close  or  near  field  condition.  In  the  near 
field  It  Is  questionable  whether  the  sound  centroid  and  the  target 
Itself  will  coincide. 

Definitive  studies  of  this  problem  have  not  been  reported  as 

yet. 

The  following  sections  of  this  report  will  describe  computer 
programs  that  simulate  the  real  time  processor  (RTF)  and  the  fast 
Fourier  transform  processor  (FFT)  algorithms  and  presents  results  of 
simulation  runs  for  various  noise  and  signal  situations. 
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3.0  MODELS  USED  FOR  SIMULATION  OF  THE  PROCESSORS 

3.1  The  RTF  Model 

The  system  diagram  used  to  model  the  real  time  processor  algorithm 
Is  Illustrated  In  figure  2.  The  Incoming  acoustic  wavefront  Is  de- 
tected by  the  microphone  array  and  Is  filtered  by  the  80  to  120  Hertz 
filter  to  eliminate  as  much  background  noise  as  possible.  This  filter 
design  was  based  upon  the  acoustic  signature  for  several  tanks  which 
were  monitored  by  the  Redstone  facility  during  the  summer  of  1977. 
Experimental  data  taken  during  that  period  was  correlated  with  acoustic 
signature  data  reported  by  other  agencies  and  a center  frequency  of 
100  hertz  with  a 1 20  Hertz  passband  was  deemed  to  be  fairly  optimum, 
when  one  considers  the  possible  frequencies  of  the  source  and  the 
possible  variation  of  those  frequencies. 

Figure  3 Illustrates  the  simulated  bandpass  filter  response 
which  has  a time  response  of 

h(t)  - 2(BW)  jslnc^CBWXt- t^)|Jjcos[(0Q(t- to)jj 


where  BW  ■ bandwidth  of  the  bandpass  filter 


and 


sine  z ■ 


sin  ff  z 
It  z 


The  phase  difference  between  the  two  microphones  Is  measured  by 
using  a phase  locked  loop  (PLL)  circuit.  As  Input  noise  and  changer 
In  the  Input. frequencies  occur  the  phase  Information  from  the  PLL 
will  vary.  These  short  term  errors  may  be  averaged  out  to  some 
extent  by  using  an  algorithm  for  short  term  averaging.  Such  an 
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algorithm  has  been  proposed  by  the  author  and  is  detailed  in  the 
following  discussion. 

Signals  in  the  frequency  range  of  100  Hertz  are  low  enough  in 
frequency  to  allow  the  phase  lock  loop  to  track  minor  i 20  % variations 
in  base  frequency  and  phase  in  a time  period  negligible  compared  to 
the  period  of  the  100  Hertz  frequency  (microseconds  versus  tens  of 
milliseconds) . A short  term  average  may  be  computed  by  storing  the 
frequency  and  phase  difference  of  the  incoming  waveforms  onlce  a 
period  for  four  or  five  periods  (approximately  40-50  ms.)  and  com- 
puting the  average  of  these  frequencies  and  phase  differences. 

This  procedure  will  average  out  the  short  term  variations  of  the 
input  waveform  to  some  extent.  How  well  this  procedure  works  may  be 
Judged  by  the  results  for  the  RTF  which  compare  very  favorably  with 
the  FFT  processor.  These  results  are  described  in  sections  4.0  and 
6.0. 

Simulation  of  the  data  Inputs  from  the  microphones  is  accomplished 
by  the  use  of  a 1 volt  peak  cosine  wave  of  100  Hertz  base  frequency. 
Noise  is  added  to  this  data  input  by  use  of  a random  number  generator 
which  is  used  to  generate  random  amplitudes  varying  from  +1  to  -1. 

These  numbers  are  generated  on  a uniform  distribution  and  are  used  as 
additive  noise.  That  is  when  the  signal  is  generated  and  sampled  for 
use  in  the  program  the  noise,  which  has  been  generated  is  discrete 
steps,  is  added  to  the  signal  samples  thereby  producing  a corrupted 
signal.  Figures  4 through  7 illustrate  the  signal  plus  noise  for 
various  noise  levels  before  and  after  filtering.  The  simulated 
signal  plus  noise  may  be  illustrated  mathematically  by: 
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f(t  + A(|>)  = 1 sin  2ir((100  + Af)t  + A(|)]  + B . 

Signal  frequency  variation  has  been  modeled  by  using  a dither 
frequency  term,  Af,  added  to  the  main  signal  frequency  f^.  Thus 
fQ  + Af  is  the  signal  frequency  used  and  Af  may  be  made  to  vary  at 
any  desired  frequency  range. 

B is  the  random  noise  term  and  is  uniformly  distributed  between 
levels  b and  a . The  rms  value  of  such  a distribution  is  given  by^ 

a - 

2/T 

when  a and  b are  ‘he  upper  and  lower  limits  of  the  noise  variations. 

If  the  signal  is  a 1 volt  peak  sine  wave  then  its  rms  value  is 
.707  . Correspondingly  if  the  noise  is  uniformly  distributed  between 
i 5 volts  then  its  rms  value  is  2.886  and  the  signal  to  noise  ratio  is 

■ rife-  ■ : .245  . 

For  the  simulation  runs  conducted  a frequency  variation  of  10.24 
Hertz  (100  Hz  to  110.24  Hz)  at  a 100  Hz/sec  rate  was  used  for  the  0.1 
second  averaging  RTF  system. 

The  output  g(t),  figures  5 and  7 may  be  found  by  convolution  of 
f(t)  and  h(t)  through  multiplication  of  their  fourier  transforms  and 
Inverse  transforming  the  result  to  obtain  g(t).  One  must  be  careful, 
however,  to  choose  a proper  sampling  period  and  sampling  Interval  so 
that  the  result  for  g(t)  does  in  fact  carry  the  information  needed  to 
represent  g(t)  in  real  time. 

Since  f(t)  is  composed  of  a periodic  signal  plus  noise  f(t)  is 
not  completely  periodic.  The  FFT  transform  however  forces  one  to  pick 
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a period  over  which  to  average  and  the  function  is  assumed  periodic 
with  that  period.  As  a result,  the  reproduction  of  g(t)  by  this  method 
may  not  contain  the  non-periodic  components  of  g(t). 

A better  method  of  simulating  g(t)  would  be  the  use  of  the  dis- 
crete convolution  model.  Since  the  bandpass  filter  is  physically 
realizable  and  passive,  its  fourler  transform  may  be  determined  analy- 
tically and  hence  its  Impulse  response  determined. 

By  computing  g(t)  through  the  discrete  convolutional  relationship 

t 

g(t)  - I f(K)  h(t  - K) 

K-0 

a more  faithful  reproduction  of  g(t)  will  be  obtained.  The  disadvantage 
of  this  method  of  simulation  is  the  additional  computation  required. 

For  Instance,  it  has  been  shown  (Stockham,  T.  G. , "High  Speed  Convo- 
lution and  Correlation,"  1966  Joint  Computer  Conference)  that  for  28 
samples  the  FFT  method  is  faster  and  for  4096  camples  the  FFT  method 
is  80  times  faster.  But  with  non-periodic  Inputs,  the  FFT  method  does 
not  always  produce  accurate  answers. 

In  programming  the  discrete  convolution  equation  a method  was 
devised  that  has  good  speed  and  a mlnlmi’m  of  memory  required.  If  one 
Is  not  careful  the  memory  requirements  can  be  as  high  as  M x N,  where 
Mis  Is  the  number  of  samples  of  the  filter  function  and  N the  number 
of  samples  of  the  signal  plus  noise  function  used.  For  2048  samples 
of  each  the  memory  requirements  would  be  greater  than  4,194,304  words 
plus  associated  bookkeeping  memory  requirements.  Obviously  this  is 
rather  high  (too  high  for  most  machines  without  resorting  to  time  con- 
suming tape  storage) . A method  was  devised  that  has  essentially  the 
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1 


same  computational  speed  and  requires  only  M + N (4096  words  for  2048 

I 

samples)  plus  associated  bookkeeping  memory  requirements. 

The  PPL  phase  comparator  which  will  be  modeled  Is  a type  II  digital 
memory  edge  controlled  network  phase  comparator.  This  type  of  phase 
comparator  In  this  application  will  stay  In  lock  over  the  full  range 
of  the  (essentially  from  5 Hz  to  5 KHz) . Phase  comparison  Is  made  i 

by  cross  referencing  the  two  channels  and  using  up-down  counters  to 
derive  an  8 bit  digital  word  which  will  Indicate  the  phase  difference 
between  the  two  channels.  The  phase  Information  Is  presented  from  a 
counter  In  digitized  form  In  approximately  1.41  degree  Increments. 

Simulation  of  the  phase  lock  loop  Is  accomplished  by  calculation 
of  the  frequency  of  the  waveform  coming  from  the  filter  h(t).  This 
determination  of  the  frequency  Is  found  by  averaging  the  periods  of 
the  first  six  peaks  In  the  output  g(t)  which  Is  created  by  f(t)  as  the 
Input.  This  average  frequency  Is  used  In  the  calculation  of  the  target 
angle. 

The  average  phase  shift  Is  determined  by  comparing  the  shift  of 
the  main  four  peak  values  of  the  output  g(t  + A41),  which  Is  created  by 
the  signal  Input  f(t  A^)  as  the  Input,  relative  to  the  main  four 
peaks  of  g(t)  created  by  f(t).  This  average  phase  shift  Is  also  used 
In  calculating  the  target  angle. 

3.2  The  FFT  Model 

The  system  diagram  used  to  model  the  fast  Fourier  transform  pro- 
cess Is  Illustrated  In  figure  8.  The  filter  bandwidth  of  50-250  Hz  was 
chosen  so  as  to  provide  a frequency  window  In  the  frequency  domain  of 
the  signal  desired.  This  bandwidth  could  be  restricted  to  a more  narrow 
range;  however,  the  system  we  are  Interested  In  Investigating  utilizes 
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50-250  Hz  bandwidth  and  hence  the  range  was  chosen  for  simu- 
lations. 

The  frequency  domain  plot  will  contain  discrete  frequency  points 
(sometimes  the  nomencloture  'frequency  bias'  Is  used)  and  the  spacing 
of  these  points  (the  width  of  the  frequency  bin)  is  Inversely  pro- 
portional to  the  period  used  for  sampling.  If  a 2 second  sample 
period  Is  used  the  frequency  points  are  0.5  Hz  apart. 

The  number  of  samples,  N,  taken  per  period,  T,  determines  the 
sampling  Interval,  t = T/N,  which  in  turn  relates  to  the  width  of  the 
frequency  spectrum  calculated.  This  relation  is  NT  and  for  the  two- 
second  sampling  period  of  f(t)  and  512  samples  per  period,  the  fast 
fourler  transform  will  yield  512  frequency  terms  with  each  term  .5  Hz 
wide.  Thus,  a frequency  spectrum  plot  would  cover  from  0 to  256  Hz. 

(In  a FFT  half  of  the  terms  generated  are  complex  conjugates  of  the 
others. ) 

The  outputs  of  the  fast  fourler  transform  devices  FFTl  and  FFT2 
are  thus  generated  by  first  using  discrete  convolution  to  calculate 
the  functions  g^(t)  and  g2(t),  and  then  a discrete  FFT  calculation  Is 
used  to  generate  FFTl  and  FFT2. 

The  outputs  FFl  and  FFT2  are  complex  numbers,  one  for  each  fre- 
quency bln.  These  complex  numbers  yield  the  amplitude  of  that  frequency 
bln  and  Its  phase  angle. 

The  amplitude  Is  proportional  to  the  energy  content  of  the  signal 
at  this  frequency.  Thus  if  we  have  a strong  signal  energy  at  a specific 
frequency  such  as  100  Hertz  and  a very  low  noise  background,  we  would 
expect  to  see  a frequency  plot  that  was  very  small  except  at  100  Hertz 
where  we  would  expect  to  see  a large  amplitude  component.  Figure  9 
illustrates  a frequency  plot  for  such  a signal. 
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Figure  9.  FFT  Simulation  Output  for  100  Hertz  Signal  (S/N~  1.225) 
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Large  amplitudes  In  specific  frequency  bins  are  selected  and  the 
phase  angle  of  these  terms  are  calculated.  The  phase  shift  of  similar 
frequency  components  is  then  calculated  by  comparing  the  phase  angle 
of  a specific  frequency  bin  for  microphone  1 with  the  phase  angle  of 
the  same  specific  frequency  bln  for  microphone  2.  The  difference  in 
the  phase  angle  is  proportional  to  the  phase  shift  in  that  frequency 
term. 

This  technique  is  sufficient  to  simulate  the  system  for  these 
signal  data  Inputs  to  be  used  and  to  compare  the  same  basic  information 
from  the  FFT  and  RTF  systems.  In  general,  however,  one  would  refine 
the  Information  from  the  FFT  to  allow  target  discrimination.  This  may 
be  done  as  described. 

The  frequency  bins  are  organized  into  groups  by  using  a grouping 
algorithm  which  applies  the  principle  that  all  frequency  bins  originat- 
ing from  the  same  direction  in  space  will  have  essentially  the  same 
phase-difference  slope.  The  phase-difference  slope  is  defined  as: 

, phase  difference  ^“^i 

Phase-difference  slope  “ ^ — » — — 

frequency  fi 

where  i identifies  the  1th  frequency  bln. 

The  grouping  is  accomplished  by  identifying  which  slopes  fall 
within  a direction  window  in  space. 

The  grouping  algorithm  step  provides  two  important  sources  of 
I data.  It  identifies  the  frequency  lines  coming  from  each  target 

I detected,  which  permits  applying  judgment  as  to  the  type  of  target 

I present,  as  well  as  provides  data  to  determine  the  best  estimate  of 

the  target  direction. 
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The  best  estimate  of  target  direction  is  computed  by  averaging 
the  phase-difference  slopes  making  up  each  group. 

1 ? ^"^i 

Average  slope  per  group  = — ) — — 

" j=l 

The  target  azimuth  is  then  computed  using  the  average  slope  per  group 
and  trigometric  relationships. 


Azimuth 


0, 


arcsin 


sin  g 
cos  Og 


where 


a = 


ji 

n 


r ^ii  2 

arcsin  l for  the  microphone  pair 


k k 


The  subscript  j refers  to  the  target  (group)  number 

The  subscript  i refers  to  the  frequency  bin  number 

speed  of  sound  in  feet  per  second  divided  by 
the  microphone  separation  in  feet 

frequency  in  Hertz  for  frequency  bin  i and 
target  j 

number  of  phase-difference  slopes  in  the  group. 


Target  discrimination  is  then  provided  both  by  the  approach  used  and 
by  performing  additional  checks. 
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4.0  RTP  AND  FFT  SIMULATION  RESULTS  ! 

i 

The  simulation  of  the  RTP  provides  an  insight  into  the  possible  | 

ij 

performance  of  such  a system.  j 

i 

Table  1 Illustrates  the  performance  of  the  RTP  for  various  con-  | 

stralnts  within  the  simulation  programs.  For  this  table  the  random  | 

i 

numbers  used  for  the  noise  amplitude  and  for  the  algebraic  sign  were  j 

drawn  from  two  separate  random  number  generators.  If  different  random 
number  generators  are  used  for  the  program  while  everything  else  Is 
held  constant,  the  Individual  results  will  vary  but  the  average  error 
should  remain  essentially  constant.  That  Is  to  say  for  one  set  of 
random  numbers  generating  noise  the  angular  results  (for  40.0°  actual 
target  angle)  may  compute  to  be  41.1654°,  while  for  a second  set  of 
random  numbers  the  computed  angle  may  be  38.710°.  The  results  vary  on 
an  Individual  basis  about  the  nominal  40.0°  but  the  errors  are  within 
.3115%  of  each  other. 

The  data  In  Table  1 Illustrates  several  points.  First  of  all, 
the  fastest  sample  rate  (.00005  seconds  per  sample)  results  In  target 
angles  that  are  generally  closer  to  the  actual  target  angle  than  the 
slower  sample  rates;  this  would  be  expected.  Also  true  Is  the  fact 
that  for  good  results  at  small  target  angles,  the  faster  sample  rates 
are  necessary.  In  fact,  the  quantization  error  due  to  slow  sampling 
becomes  prohibitive  at  .00025  seconds  per  sample  time. 

For  the  simulation  areas  comparing  the  FFT  and  RTP  systems,  the 
sample  rate  used  was  .00005  seconds  per  sample  time. 

In  general,  the  results  are  Indicative  of  a system  algorithm  that 
would  be  suitable  for  a tracker,  especially  for  Initial  lockon  and 
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tracking.  As  the  vehicle  nears  the  target,  an  IR  tracker  with  a 
narrow  field  of  view  could  be  used  for  the  terminal  homing.  This 
would  most  likely  be  advisable  since  the  IR  tracker  could  provide  a 
faster  control  loop  in  the  guidance  system. 

From  the  data  of  Tables  1 and  2,  it  is  evident  that  small  varia- 
tions in  the  phase  shift  and  the  inverse  of  frequency  result  in  signifi- 
cant errors  for  small  target  angles.  This  is  to  be  expected  since  the 
angle  calculation  is  proportional  to  the  Inverse  cosine  of  the  phase 
shift  and  the  inverse  of  the  frequency.  For  small  angles  the  cosine 
function  has  the  fastest  rate  of  change  as  opposed  to  the  larger  target 
angles  which  correspond  to  the  slowest  rate  of  change  of  the  cosine 
function  versus  the  argument  of  the  function.  Another  way  of  saying 
this  is  to  state  the  change  of  the  argument  of  the  Inverse  cosine 
function  is  slowest  around  the  90°  area  (the  broad  peak  of  the  function) 
and  quickest  around  the  0°  area  (the  steepest  slope),  which  may  be 
Interpreted  as  implying  that  small  changes  in  the  argument  of  the 
inverse  cosine  function  yield  proportionally  larger  changes  in  the 
resulting  angle  calculations  for  arguments  corresponding  to  very  small 
angles. 

In  any  event,  the  results  indicate  (as  we  know  to  be  true)  that 
both  types  of  algorithms  make  suitable  null  trackeLS  by  keeping  the 
target  broadside  (6  * 90°)  to  the  axis  of  the  array  which  operates  in 
the  most  accurate  region  of  the  array. 
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5.0  COMPARISON  OF  THE  TWO  PROCESSORS 

The  comparative  results  of  the  performance  of  the  RTP  and  FFT 
processor  algorithms  are  depicted  In  Table  3.  The  data  In  this  table 
was  made  using  a different  random  number  generation  method  than  that 
used  to  generate  the  data  of  Tables  1 and  2.  However,  the  same  random 
number  generator  was  used  for  both  the  RTP  and  the  FFT  simulations  of 
Table  3. 

In  this  comparative  run  the  two  algorithms  were  used  as  they  would 
be  In  an  actual  system.  That  Is,  the  RTP  system  was  presented  with  the 
signal  plus  noise  data  and  an  answer  was  computed  In  . 15  seconds 
(actually  the  answer  Is  computed  using  .08  seconds  of  the  data  emanat- 
ing from  the  filter  starting  at  t = 0),  as  opposed  to  the  two  seconds 
required  by  the  FFT  system. 

The  same  signal  data  and  the  same  noise  data  were  used  for  both 
the  RTP  and  the  FFT  calculations  and  the  results  Indicate  the  relative 
performance  under  essentially  the  same  Input  conditions. 

Good  agreement  Is  obtained  for  both  processor  algorithms  although 
the  FFT  processor  Is  the  more  accurate  of  the  two.  In  fact,  for  large 
noise  levels,  the  FFT  performs  much  better  than  the  RTP,  which  was 
unable  to  arrive  at  a satisfactory  answer  for  small  target  angles. 
However,  If  both  systems  are  around  the  null  point,  both  are  useable 
for  tracking.  Another  point  to  consider  Is  the  fact  that  as  the 
vehicle  approaches  the  target  the  signal- to-nolse  ratio  will  Increase 
and  likewise  the  accuracy  of  the  processor  algorithm  Increases  In  both 
cases. 

Table  4 Illustrates  the  results  of  a simulation  In  which  the 
frequency  of  the  Incoming  signal  Is  varying.  For  this  situation  the 
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FFT  processor  will  suffer  a degradation  In  performance  that  would 
render  It  unusable  for  a null  tracker  unless  a considerable  addition 
of  software  processing  Is  added  to  the  algorithm. 

Figure  10  Illustrates  the  spectrum  which  the  FFT  system  "thinks" 

It  sees  when  presented  with  a varying  signal  frequency.  Since  the  FFT 

uses  all  data  over  the  2.0  second  Interval,  It  "thinks"  that  energy  ! 

exists  over  the  range  of  frequencies  that  the  signal  traverses  during 

this  time. 

By  comparison  the  RTF  system  holds  up  well  since  the  signal  fre- 
quency does  not  vary  appreciably  In  the  short  processor  time  required.  ' 

For  this  simulation,  the  signal  frequency  was  varied  at  a 100 
Hertz  per  second  rate  for  the  RTF  system  and  27.90  Hertz  per  second 
rate  for  the  FFT  system. 

Finally,  In  Table  5,  are  simulation  results  of  the  FFT  system  for 
a 0.5  second  averaging  time;  which  was  proposed  for  a quicker  updating 
system  for  close- In  tracking. 

Results  Indicate  degradation  of  the  system  performance  with  short 
sample  time  for  a constant  signal  frequency  the  results  are  useable 
for  a high  slgnal-to-nolse  ratio. 

The  results  for  a single  frequency  that  varies,  however,  Indl- 
.:ate  that  the  system  would  not  be  useable  since  the  output  results  are 
very  erratic  In  accuracy  (In  fact,  almost  constant  angle  output 
regardless  of  target  angle) . 
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6.0  CONCLUSIONS 

Computer  programs  have  been  written  that  simulate  real  time  pro- 
cessing (RTP)  and  frequency  domain  (FFT)  algorithms  for  processing 
acoustic  data  such  as  that  emanating  from  tanks.  Simulation  tests  of 
these  systems  have  been  accomplished  and  the  results  summarized  In 
this  report. 

The  simulation  results  Indicate  the  following  points: 

1.  The  FFT  system  Is  the  most  accurate  under  conditions  for 
which  the  signal  frequency  Is  constant. 

2.  The  RTP  system  will  provide  a good  tracking  signal,  especially 
for  target  angles  of  30®  or  more. 

3.  The  FFT  system  cannot  tolerate  a signal  with  varying  frequency. 

A.  The  RTP  system  Is  capable  of  processing  the  data  approximately 

10  times  faster  than  the  FFT  system. 

It  seems  feasible  to  design  an  acoustic  tracking  algorithm  that  uf^lzes 
a real-time  processing  algorithm.  An  RTP  system  holds  several  advantages 
over  an  FFT  system,  namely: 

1.  It  Is  much  less  expensive  due  primarily  to  the  simple  hard- 
ware required  with  respect  to  an  FFT  system, 

2.  It  Is  an  order  of  magnitude  faster  In  processing  time  that 
the  FFT  system  which  yields  a faster  guidance  system  control  loop. 

3.  It  has  the  ability  to  process  and  output  useable  data  for 
signals  with  signal  frequencies  that  are  varying. 

It  Is  recommended  that  a real-time  processing  algorithm  for 
processing  acoustic  signals  such  as  those  emanating  from  a tank  be 
prototyped  In  breadboard  form  for  evaluation  In  an  actual  environment 


under  laboratory  and  controlled  field  test  conditions. 

It  Is  further  recommended  that  software  be  developed  that  Is 
suitable  for  a microprocessor  to  average,  categorize  and  perhaps 
provide  a "smart"  tracking  algorithm  for  the  RTF.  Using  this  soft- 
ware development  or  tracker  system  using  the  RTF  approach  may  be 
designed  for  a prototype  tracker  breadboard  development. 
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APPENDIX  A 

The  following  Is  a computer  program  written  in  fortram  for 
simulation  of  the  real  time  processor  algorith,  RTP.  This  program 
was  run  on  a UNIVAC  1107  mod  II  system. 


i 


i 


I 


THIS  PKOGRANi  IS  THL  SOLUTION  TO  THE  REAL  TIME  PROCESSOR 

H(T)  IS  the  impulse  reponse  of  the, system 

XlT)  IS  THE  INPUT  SIGNAL  TIME  FUNCTION 

these  parameters  Are  to  be  set  for  fach  case 

pTbAE  = pekCent  Target  angle  error 

sFo  = SIGNAL  Frequency 

SPO  = SIGNAL  PERIOU  = 1/FKEOuFNCY 

PnSF  = actual  PHASE  SHIFT 

AlARG  = actual  TARGET  AnGLe 

ut»  F»  G»  IFQ 

eixter  ID  the  filter  bandwidth  in  hertz 

ENTER  IFO  the  FILTLR  CEnTEr  FREQ  IN  HErTZ 
EhTER  UT 

G must  always  be  LESS  THAN  OR  EQUAL  TO  F 

enter  g the  « of^ut  increments  in  H(T) 

enter  F»  the  t*  OF  UT  increments  IN  xU) 

enter  the  MiOPTS  AmPTS  OF  OT  FOR  EACH  DT  IN  H(T) 

enter  the  MIDPT  AMPTS  of  UT  FOR  EACH  DT  IN  XT) 

Y»V  MUST  BE  dimensioned  TO  (Q+F) 
integer  g»f 
integer  Z 
integer  R 

G=1024 
F=i;04B 
DT=. 00005 


40 


L 

c 

y. 

t 


110 


oOb 

o04 


oo;^ 

b03 

uOb 


bOO 

34 


bOl 

35 

ill 

C 

c 


^b 


13 


OlMENblOr  P^YWM(400)  tPKYWUOu)  » PKYM  ( 400 ) »PK  Y ( 400  ) 
OlMENblON  XK(35fl4) »X(35a4) » Y»(4096) 

DlMENblON  hl2048) »v(35a4) 

VLb  = VELOCITY  OF  bOUNU  IN  FlET  PER  SECOND 
VLb  = 1100. C 

oaky  = distance  BlTWEEN  MICROPHONES  ON  ARRAY  IN  FEET 
oaky  = 4.0 
SFO=lCiO.O 
SPo  = 1/SFO 

THIS  SETS  NOibL  To  0 IF  INb=0  OR  GIVES  NOISE  IN  INS=1 
INb-O • 0 

PhsF=123. 01430 

THIS  oENERAIES  the  DT  MIDPT  VALUES  OF  X(T)  AND  H(T) 

PI  = 3.1415920 
JK  = G+  F 
Lll  = G+F 
lu=40^ 

ifo=igo 

DO  llO  JV=1»G 
VluV)=OV 

VV=(2»(JV-(0/2) j+l)*0T/2f 

ri<JvT=(2*IL)*SlN(Pl*lB*Vj)»C0S<2*PI*lF0*VJ)  )/(Pl*lB*VJ) 
continue 

Call  PL0TS2‘V*H»G»1) 

00  117  IJV=1.2 

- 1.0 

0 
0 
0 


) GO  TO  602 


10  = C 
Jt  = 1, 

10  = 1, 

11  = 1, 

JvJmO  * 0 

IKINb  .GE.  1.0)  00  TO  604 
AiNS=lOOOOo*0 
00  60b  1=1, F 
XK( I )=0.0 

continue 

GO  TO  606 

continue 

A1NS=1 .0/lNb 
IrdJV  .EO.  1 
XK( 1 )=385. 

GO  TO  603 
XK{1)=I30. 
continue 
call  hanuu(ar 
continue 

DO  111  KV=1»F 
VK=(2*  (KV-l)+l)*p..  ... 

PHI  = (2*PI’>PHSFT/360. 

IF(IJV  .GT.  1.)  GO  TO  33 
00  TO  34 

IKXR<F-KV)  .LT.  0,5)  GO  To  oOO 
X(KV)  = lCOS(i:*PI*(bFO+KV*.005)*V 

^?=icoSt,i*PI*(bFO+KV*,005)*VK)  )-(INS*XR(KV)  ) 
5^ 


pF) 


)r./2. 


VK) )+( INS*XR(KV) ) 

xTkvT:"". 

GO  TO  35 

iF(XRtF-KV)  .LT.  0. 

XlKV)=(C0S( T2*Pl*Tb 
GO  TO  35  , 

XlNV)  = (C0Sl  (2*Pl*lbFO-»KV*.005)*VK)+PHI)  )-(  INS*XR(KV)  ) 
continue 
continue 

INSERI  the  Do  LOOP  AND  PLOTTING 
This  reverses  the  ORDER  OF  XTT) 

R=F/2.0 
00  25  1=1, R 
T=X( I) 

Xll)=X(Z) 

XtZ)=T 


J.5)  GO  To  601 
ISFO+KV*.005)*VK)+PHI) )♦( INS*XR(KV) ) 


ROUTINE  FOR  X HERE 


Continue 


i9-!$  = 
y=o.o 


= 1»JK 


DO  . 

Y«(K.  - - 
IKIJ  .LF 
GO  TO  14 

continue 

DO  20  I=l»Io 
K1=(F*1.0)-1 
Kt=(IJ>1.0)-I 
tOM=H(KY)*X(Kl) 


O)  60  Tu  13 


♦ DT 


SRX8  PAQI 18  BBST  dUALITY  PBACIICABLE 
ffpM  nM>Y  gUBMlSHSD  rODD^g  — ^ 


41 


tiO 

14 


15 

16 


tl 

17 

10 

C 

525 

C 


Jl4 

513 

C 


29 

c8 


421 

126 

C 

c 

117 

319 

j2U 

321 

322 
C 


Y*»(K)=£:UM+Y»*(i<) 

CUNTHHJE 
lo=IJ+l .0 
bO  TO  17 

continue  , 

IKJJ  .GF.  tp-G))  GO  TO  16 

DU  15  JL=l»u 

K1=(F“JJ)-0L 

Kr=(6+1)-JL 

tUfi/izHlKYj+ACKD^DI 

Y»*  (K)=EUM-*-Y>n  (K) 

continue 
jJzJJ+1.0 
GO  TO  17 
Continue 
11=11+1.0 
K1=(F+1)-I1 
Do  21  1=1 » lu 
KG=G+1.0-I 
00=16-1+1 

t;ONi=H(KG)+X(jG)*Dl 
Y**(K)=EUV  + yw(K> 

continue 

I0=IG-1 

continue 

coi,tinue 

PUT  PLOT  KOUTINE  PoP  Yw  HEHL 
DO  325  JIV=1»LLL 
V(oIV)=JlV 

continue 

C«LL  HL0TS2tV» YW»LLL»li 
This  PIND5  IHL  PEakS  OF  YW 
Ko=l.C 

PNYWMll)=Ywll) 

PKYW<1)=1.0 
DO  313  J=2»LlL  . 

IKYWU)  .GT.  YW(o-l))  GO  TO  313 
If ( J .EG.  2)  60  TO  313, . ^ 
1+(YW(J-1)  .gT.  Y»i(J-2))  Go  to  314 
GO  TO  313 
PNYWMfKJ)=Yn(J-l) 

PAYW (K0)=J-1 
Ko=KJ+ 1 
KKo=Ko 

Continue 

This  finds  IhE  main  PEAK  OF  YW 

X=pKYrtM(l) 

MPKYW=1 
DO  2«  J=2»KKJ 
Y=PKYi»M(J)  ^ 

IP ( Y .GE.  X)  GO  To 

GO  TO  2fl 

X=Y 

MPkYW=J 
COrjTINUF. 

IPIIJV  .gT  _ 

DO  421  J=1»NKJ 
PKY ( J)=PKYW Ij)  . 
pNYM( J)=PKYWMl J) 

Continue 

mpky=npkyw 

JoK=KKJ 

continue 

insert  DO  LOOP  and  WRITE  ROUTINES 
Y IS  IffE  PHASE  shifted  SIGNAL  AND 
COtjTiNUE 
WhlTE(6»319) 


1.0»  60  TO  126 


FohMAI ( IHl » ‘PEAK^VALUES 


»KltE  ( 67320 )'  TpKYliT jy7;^PKYwM?jf?J=i  »KKJ) 
-OhMAT (4(G12,7»2X»G12«7) ) 


F(T)») 


-ott^TtlHoJ’PEAK  values  OF  Y EQUAL  TO 
•KITE 16» 322)  TPKY(d) fPKYMiJ) *3=1*00X1 


F(T+PH1) •) 


FohMAT (4T612.7*2X*612.7) 1 
SDp=0*0 

THIS  PiNDS  Average  period 

DO  317  0=1*4 
K6=MPKYW-2+0 
0«=MPKY«l-3+o 
DlFF  = PKYwfKU)“PXYW(00) 


USING  2 PEAKS  EACH  SIDE  OF  MAIN  PEAK 

nzi  ?AQI  ZS  Basi  quality  FRACnCABLt 
nmi  ooPY  ruRHisHSD  xolqq 


1 


42 


J17 


J16 

C 

bO 

61 


o20 

b2 


o3 

44 

45 
C 


40 


43 

•♦2 

J2b 

327 

323 

324 


SOM=SUM+DIFf-' 

continue 

AvHRD=(SUM/4.0)*m 
FKlQ  = 1/A^PrU 
bMFAS=0.0 

THIS  unds  Average  phase  shift  using  2 peaks  each  side  of  main 

DO  3lb  J=l»4 

K0=MPKy-2tj 

J0=MPKYW-2+J 

prtS0F=(  (PKYtKO)-PKYW(  JQ)  ) *i..T*360 . 0 ) /AVPRD 
SWPAS=SVPAS  4 PASDF 

continue 

THIS  CHECKS  THE  AVPAS  FqR  360  ROTATION 
IHSFHAS)  bU, 61*61 
AvpAS=ABS(SHpAS)/4.0 
GO  TO  62 

IFISMPAS  .LL,  1440,0)  Go  To  b20 
SPPAS=SWPAS“1 440.0 
GO  TO  61 

AVPAS=360.U'(SMPAS/4.0) 

CONTINUE 

IFIAVFAS  .LE,  360*0)  60  TO  63 
AwpAS  = AVpAs  - 360.0 
GO  TO  62 

continue 

wRITE<6.44) . 

format  (1110.  'mPKYW*  . 12X.  *MPkY*  . 13X.  * AVPrd*  » l?X*  • AVPAS* ) 
WHITE(6.45). (MPKYh.MPKY.AVPRD.AVPAS) 
format  (G12. /»bX*Gli;.7*5x*Gl2.7,5X*Gl2.7) 

THIS  PROTECTS  AGAINST  INVALID  ACOS  ARGUMENTS 
ARt>VT=(AVPAS*VLS*MVPRO)/(DARY*360.0) 

UV=ARb(ARGMT ) 

IF (RV  .LT.  1.0)  GO  TO  40 

ANoLE=9g9.99 

AlARG=Q99.y9 

Plt.AE=9g9.99 

GO  TO  42 

ANgLE=(ACOS(BV) )*(360.0/(2*P1) ) 
ATaMT=(PHSF»VLS*SPo)/(DaRY*360.0) 

Bd=AfiS(ATAMT ) 

IF  IBR  .LT.  1.0)  GO  TO  43 
ATARG=999.y9 
PI AGE=999.y9 
GO  TO  42 

).0/(2*pi)  ) 


AIaRG=(ACOS(ATAMT) )*(360.0/  , 

PIv3AE=(  (ATARG-ANGLE)/ATAR6>*100.0 


continue 

wRITE(6*32o) 

Format (iHO* ’RTp  with  s/n  = i/Ins*) 
wH1TE(6*327)  AINS 
format (G13.0) 

WR1TE(6*323) 


frequency '.SX. ’AVG 


. SHIFT* »5X* .actual  PHASE  SHlF 
iRGET  ANGLE* »5X»*«  TARGET  ANGL 


format  (IHO*  ' ... 

/T*  .5X*  * AMGLt.  TO  TARGET*  »5X»  .ACTUAL 
/E  error*  */) 

WRITE(6.324)  (FREU*AyPAS*PHSF*ANGLE»ATAR6»PTGAE) 

Format (e13.o,4X *Gl2. 7* 6X* 612,7* 11X»G12.7»8X* G12. 7* 13X» 612 •?) 

STOP 

ENL 
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APPENDIX  B 

The  following  is  a computer  program  written  in  fortram  for 
simulation  of  the  frequency  domain  processor  algorlth,  FFT.  This 
program  was  run  on  a UNIVAC  1107  mod  II  system. 


The  NAIwE  PKYkVM  ahpeaks  in  a OiMinSION  or  type  statement  but  Is  never  rfferfn 

THE  NAPE  PKY*I  ApPtARb  IN  A DIMENSION  OR  TYPE  STATEMENT  RUT  IS  NEVER  REFERENC 
This  PROGRAM  IS  ThL  SOLUTION  To  THE  FFT  PROCESSOR 
HIT)  IS  THE  IMPULSE  HEPqNSE  OF  THE  SYSTEM 
X(l)  IS  THE  INPUT  SIGNAL  TIME  FUNCTION 

these  parameters  are  to  be  Set  for  FACh  case 
p)gae  = percent  target  angle  Error 

SEo  = SIGNAL  FREQUENCY 

SPO  = SIGNAL  period  = I/FReQUENCY 

PhsF  = actual  phase  SHIFT 

aiarg  = actual  target  angle 

DU  F»  Gf  10,  IFO 

enter  IR  the  filter  bandwidth  IN  HERTZ 

enter  IFo  the  filter  center  FREQ  IN  HERTZ 


xaxs  FAdE  18  BEST  QUALlTTf  PBAOKOAM* 
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C 

C 

C 

C 

c 

c 

c 


c 

c 


c 

c 

c 


ilO 


QOt> 

o04 


b02 

u05 

oOo 


ENTER  UT 

G KUSl  ALWAYS  BE  LtSS  THAN  OK  EQUAL  TO  F 
ENTER  G THE  tt  OF  bT  INCRFMENTS  IN  H(T) 
enter  Fr  THE  tt  OF  UT  INCREMENTS  IN  X(T) 

El'.TER  THE  MiORTS  AMPTS  OF  UT  FOR  EACH  OT  IN  H(T) 
enter  the  MIDPT  AMhTS  of  DT,F0K  EACH  DT  IN  XT) 

Y»V  MUST  BE  dimensioned  TO  (g+F) 

IimTEGER  G»F 
IimTEGEH  2 
IiNiTEGtR  R 

IE  = THE  power  of  2 FOR  THE  FFT  ROUTINE  NUFFT(IE) 

IL=10 

X«  = THE  number  OF  NOISE  TeRmS  TL  BE  ADDED  TO  THE  SIGNAL  X 

Ih  = number  OF  points  For  FFt  to  operate  ON 

Ih=1024 

G=bl2 

F=3504 

Dr=. 000125  , , „ . 

DIMENSION  PnyWM(400)  rPKYWUOo) 
dimension  FPkFT(400) »FPKNB(400) 

DIMENSION  - ' . 

DIMENSION 
DIMENSION 
DIMENSION 
DIMENSION 


rtrv'  \ -ry  V / 

XK(3ba4) , YW(4096) »X(3584) »V(3584) 
HIiQ24) »Y(1024) 

AYwll024) ,FIY(1024) »PKFFT(40o) »PKI 
Fab(1024) ,FI (1024) 


C0S(2*PI*IF0*VJ) )/(Pl*lB*VJ) 


.PKFNR(400) 

-1024^ 

u*Mr.iNJiur4  KY\|^(1024) 

COmMON/HOOK/xREAL(1024) »XIMA6(1024) 

Ev>UlV«EENCE  ( XREAL  ( I ) » AYW  < 1 ) ) » ( X IMAG  ( 1 ) » F I ( ) ) ) 

Vi-S  = VELOCITY  OF  SOUND  IN  FlET  PER  SECOND 
VL-S  = llOO.O 

U«RY  = distance  between  MICROPHONES  ON  ARRAY  IN  FEET 
DArY  =4.0 
SFo=lU0.0 
SPu  = 1/SFo 

THIS  SETS  NUISE  TO  ZERO  iF  INS  = 0 OR  GIVES  NOISE  IF  INS  = 1*0 

INS=1*0 

PhsF=130.000 

THIS  generates  THL  DT  MIDPT  values  of  xd)  AND  H(T) 

Pi  = 3.1415920 
JK  = G+  F 
LLL  = G+F 
IFO=150 

Iti=200 

DO  lie  JV=1'G 
V(OV)=JV 

Vo=(2*( JV-(G/2) )+1)*DT/2 
Hi  JV)  = (2+IU*SlN(Pl*IB*Vj)* 

continue 

call  F^L0TS2(V»H*G»1) 

DO  117  IJV=1,2 
DO  53  J=1»1H 
FI  I J)=0.0 

continue 

IG  r b - 1.0 
JY  = 1.0 
11  = 1.0 
U = 1.0 
JorO.O 

IF(INS  .GE.  1.0)  GO  TO  604 

ains=iooooo.o 

DO  605  I=1»F 
XM( I )=0.0 

continue 

GO  TO  60b 

continue 

A1NS=1. 0/INS 

IF ( IJV  .EO.  1. ) GO  TO  602 
XK( l)=385. 

GO  TO  603 
XK(1)=130. 
continue 
call  RANDU(ar»F) 
continue 

DO  111  KV=1;F 
Vn=(2* (KV-i)+l)*pT/2. 
pHi  = (2*PI*PHSFT/360. 
l+dJV  .GT.  1.)  GO  TO  33 
GO  TO  34 

IF(XRIF-KV)  .LT.  0.5)  GO  JO  bOO 

X (KV)  = (COS(.:  + Pl*(bFO+KV*. 01395)  *VK)  ) + (lNS*XR(KV)  ) 

GO  TO  35 


XBX8  PAia  IS  BBST  QUAIilH  PBACXIOAW* 


45 


uOO 


oOl 

J5 

111 

C 

C 


13 


<20 

14 


15 

16 


<;! 

17 

10 

c 

c 

c 


03 
12b 

04 
C 

52 

C 


1»  JK 

. t.)  GO  To  13 


GO  TO  16 


X (kV  ) = (COS  (<:*HI*(bFO+KV*.  01395  )*VK)  ) - ( INS*XP  (KV)  ) 

GO  TO  35 

1F(XR(F-KV)  .LT.  Q.5>  go  TO  601 

X(KV)={COS( ( 2*Pl*(bF0+KV*« 01395 )♦ VK)+PHI) )+(INS*XR(KV) ) 

GO  TO  35 

X(kV)=(C0S( ( 2*PI*(SF0+Kv*» 01395 ) ♦VK ) +PhI ) )-(INS*XR(KV) ) 

COnTIIjUE 
COMIIjUF 

II<ibERl  THE  00  LOOP  AND  PLOTTING  ROUTINE  FOR  X HERF 
TnlS  KEVERSos  THE  ORDER  OF  X(T) 

<'=T 

R=f-/2.0 
DO  25  1 = 1. R 
T=X(I) 

X(l)=X(Z) 

X(Z)=T 
Z=<i-1 
CONTINUE 
DO  10  K = 

YI(.(K)=0.0 
IF(IJ  .LE 
GO  TO  14 

continue 

DO  20  1 = 1. lU 
K1=(F+1.0)-1 

KT=( 1041.0)71 
E0k=H(KY)*X(k1)*D1 
Yw  (K)=EUM+Yiiv  (K) 

continue 
lo=IJ+l  ,0 
GO  TO  17 

continue 

IF(JJ  .GE»  (F-G)) 

DO  15  JL=1.G 
Ki=(F-JJ)-OL 
Kt  = (G-*l)-JL 
EUm=H(KY)*a(kI)»DT 
Yw(K)=EU^^  + Yw(K) 

continue 

Uo» jU4 1 • 0 
GO  TO  17 

continue 
11=11+1. c 

Ki=(F+1)-I1 
DO  21  1=1. lo 
KO=G+1.0-I 
JO=IG“I+l 

EUN=H(KG)*XIjG)*DT 
YW(K)=EUy+Yw(K) 

continue 

IG=IG-1 

continue 
continue 

insert  WRITu  FOR  OUTPUT  YW  HlRE 
insert  do  loop  and  PLOTS  ROUTINE  FOR  Yw  HERE 
THIS  LO  LOOP  TAKES  THE  PROPER  NUMBER  OF  POINTS 
UV  63  J=1.IH 
V( j)=o 
K=4*U 

ATv,(  J)=YW(K) 

continue 

1+ (IJV  .gT.  l.O)  GO  TO  12b 
GO  TO  127 

continue  , 
call  NUFFT(lE) 

Do  64  MM=1.1h 

FYkkTv^)=SQRT  (AYW(MM)**2+FI  (Mm)**2)*DT*2 

continue 

RLaL  part  in  AYW  IMAG  part  in  FI  MAG  PART  IN  FYW 

IG=IH/2 

DO  52  0=1. lo 

V( j)=0 

ConTII^UE 

call  PL0T52  V*FYW.1G»1)  ^ 

THIS  finds  iHt  peak  VALUES  OF  FFT 

FPkF T { 1 ) =F  Y rt  ( 1 ) ]»IXS  PAO*  18  BBSI  QUALITY  PRACTICABW 

11h=  IM/2  JWQJIOOiPYBJjpttlSip)  IOI^ 

Do  6n  0=2. Hh 


FROM  YW(T)  for  FFT 
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<37 

<38 

83 

84 

il7 

C 

C 

C 


HO 


43 

42 

326 

j?7 

023 

024 


11K=MHKFT 

rtKlTE(6»87) 

FOKMAI  (1H0»  'YHIP)  * .12X,  »FiY(IIP)  •) 
W8lTE(6»e8)^(<?0l  »ZTI ) 

FOKMAT (G12.7,4X»612.7) 

PhA6L=ATAN(t I Y ( I IFT/Y ( I IP ) ) *180. 0/PI 
WPlTE (6»83) 

“ I TimO» 'mPKFT* , 12X» ‘PHaGL* ) 


wHlTE(6»fl4) 

FOKMAITg  ' 


FORMA!  

»fl4)^(MpKFT»  PHAgD 
G12.7,4X»G12.7) 

IKIJV  .gT.  1.0)  60  TO  gl9 
PhLiLY=PHAGL 

continue 

continue 

THIS  Finds  thl  fplouency 

FKtQ=‘MPKFT)/( (6+fT*DT) 

This  TindS  Tut  PHASE  SHIFT  OF  MAIN  PEAKS 
pAsDF=(PHDLt-PHAGLT 

THIS  protects  against  invalIu  acos  arguments 

AKGMT=  { PASUt  ♦ '/LS/FkEQ  ) / ( DAr  Y ♦ 360.0) 
tiW=ABS(ARGMT) 

IF (BV  .LT.  i.O)  Go  TO  40 

ANGLE=999.99 

ArARG=999.99 

PIgAE=999.99 

00  TO  42 

ANoLE=(ACOStBV) ) ♦ ( 360 » 0/ ( 2*Pl ) ) 
AlAVT=(PHSF*VLS*SPo)/(t)ARY*3toO.O) 
bU=AbS( ATAMT ) 

IF(BB  .LT.  1.0)  Go  TO  43 
AImR6=999.99 
PUGE=999.99 
GO  TO  42 

Al ARG=( ACOS ( ATAMT) >*(380.0/ (2*PI) 
pioAE=( (ATARG-ANGLl) /AT ARG) *100.0 
continue 

wRlTE(b»326) 

FORMATTihO.  FFT  with  S/N  = l/lNS’) 
WK1TE(6»327)  AINS 
forma  r 


) 


(G13.8) 
WKITElb»323) 
FORMA! (1H0» ' 
/Tssxf  * angle 
/E  lRROR*./) 
wRlTE(b.324) 

formatTg' “ 

STtP 
ENL 


FREQULNCY* »5Xi • AVG  PHASE 


[FT* »5X» 'ACTUAL  pHASE  SHJF 


TO  TARGET' »5X» 'actual  TARGEi  ANGLE' »5X» 'HCt'aRGET  ANGL 

‘FRE(>i»PASUF»PhSF» ANGLE. ATARG»PT6AE) 
G12.7,bX*Gl2.7»6X»Gi2,7.lIx»Gl2,7.8X»G12.7»l3X.G12*7) 


tni  rui  I.  BKi  auALtiY 

moil  OOf  Y lyBHlSHSC  lOWiQ 


1 


48 


SUbROLTiNE  NuFFT<KE>  ^ , 

CUn,MON/HOOK/XKEAL(1024)  » XlNiAb  ( 1024 ) 
NU=KE^ 

n=^**nu 

N*J=N/2 

nui=nu-i 

DU  100  L=1»1^U 

102  Do  101  I=l*»'‘2  , 

P=lBI  »K(K/2**tMUl»NU) 
AKo=b.2831U‘J‘*P/FL0AT(N) 

C=COSi ARG) 

S=‘oINl  ARG) 

KlzK+l^^^^ 

TKlAL=XREALJk1N2)*U+XIMAG<K1n2)*S 
TlMAG=XlMAG<KlN2)»C-XREALiKlN2)*S 
XKl  AL IK 1N2 ) =Xi^EAL  iK  1 ) “TREAL 
XiMAGUlh2)j=XIMAGtKl)-TlKAG 
XKLAL<K1)=XKeAL  K1)+TREAL 
XlKAG(Kl)=XlMAGtKl)+TIMAG 

101  k=k+i, 

K-K+N<i 

Ir(K.LT.M)  Go  TO  102 
K-0 

nui=nli-i 
100  N*i=N?/2 

DU  103  K=l»iM  ^ . 

1 = 1BITK(K71»N<J)+1 

IfIi.LE.K)  go  to  10,3 
TRtAL=XREALlK) 

Ti|wiA6=XlMAGlK)  , , , 

Xr<LAL  K)=XI<oAL(I 
XaMAGIK)=XIMaG(I ) 

XKtALlI)=TKLAL 
XlMAGl I)=TIMAG 

103  continue 
rliurn 
end 


function  ItilTR(J»NU) 


J1=J 
IblTR=0 
DO  200 


I=1»NU 


<.00 


J<r=Jl/2 

IolTR=lBlTR*2+( J1-2+J2) 

U 1 — «J2 

Kt-TURN 

end 


SBIS 


XOSD.Q 


